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Abstract. - We present a dynamic Monte Carlo study of the Kosterlitz-Thouless phase tran- 
sition for the spin- 1/2 quantum XY model in two-dimensions. The short-time dynamic scaling 
behaviour is found and the dynamical exponent 9, z and the static exponent n are determined 
at the transition temperature. 



The existence and the nature of the phase transition in the quantum XY model is a 
long-standing problem. In 1973, Kosterlitz and Thouless explained what is now called the 
Kosterlitz-Thouless (KT) phase transition in the classical XY model, in terms of topological 
order, characterized by an exponentially divergent spatial correlation length and susceptibility 
[Q. General universality arguments suggest that the same KT transition may occur in the 
quantum XY model ^, However, a quantitative determination of the critical exponents 
and the transition temperature with Monte Carlo methods is very difficult since one suffers 
from critical slowing down. 

Due to the exponential divergence of the correlation length at the transition temperature 
Tkt and the fact that the system remains critical below Tkt, numerical simulations of 
critical systems with a KT transition are more difficult than those with a second order 
phase transition. The situation is even more severe for quantum spin systems with a KT 
transition ||, || . A standard approach to the quantum XY model is the quantum Monte Carlo 
(QMC) method where the Suzuki-Trotter transformation is used to transform the quantum 
system to a classical one || [tJ. For the 2-dimensional spin- 1/2 quantum XY model, Loh, 
Scalapino and Grant first estimated the KT transition temperature between Tkt — 0.4 — 0.5 
[§. The authors of refs. |, § improved the results with extensive QMC simulations on 
lattices up to 128 2 . With a loop-cluster algorithm |], [Io| which is often more efficient than the 
conventional QMC methods, Harada and Kawashima recently measured the helicity modulus 
for temperatures between T = 0.2 — 0.60, and determined rather accurately the transition 
temperature Tkt = 0.3423(2) on the lattice 64 2 pi. In this loop-cluster algorithm a loop 
is formed by spin-pairs on the interacting plaquettes and all spins on the loop are flipped 
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simultaneously to overcome critical slowing down [ |12[ . However, up to now it is still difficult to 
determine the critical exponents in equilibrium accurately. For example, the critical exponent 
r\ has been given as rj = 0.25 ± 0.01 in ref. [|, and r\ = 0.290 ± 0.09 or rj = 0.276 ± 0.014 in 
ref. ||. These values differ from each other. 

On the other hand, in recent years much progress has been made in critical dynamics. 
Traditionally it was believed that universal scaling behaviour only exists in equilibrium or in 
the long-time regime of the dynamic evolution. However, for classical magnetic systems, it 
was discovered that universal scaling behaviour emerges already in the macroscopic short-time 
regime 0, [L5|, |l6|, [L7|, [L8|, |ig|, |(| [H], II- Important is that new independent exponents must 
be introduced to describe the dependence of the scaling behaviour on the initial conditions, 
or specify the scaling behaviour of special dynamic observables. More interestingly, based on 
the short-time dynamic scaling form, static exponents and the dynamic exponent z originally 
defined in equilibrium or in the long-time regime of the dynamic evolution can be extracted 
already from the universal short-time behaviour Jl^, ^3, 24 . This provides a possible new way 
out of critical slowing down. 

In this paper, we investigate whether there exists universal short-time dynamic scaling 
behaviour in quantum spin systems, taking the spin-1/2 quantum XY model in two dimensions 
as an example. We determine the new exponent 9 and the dynamic exponent z as well as the 
static exponent 77 from the power law behaviour of the observables at the beginning of the 
time evolution. 

The spin-1/2 quantum XY model is defined by the Hamiltonian 

H=-J E[<" s i+<-4] > (!) 

<ij> 

where < ij > stands for nearest-neighbour pairs on a 2-dimensional lattice, sf and s? are 
spin operators defined on each lattice site, which can be expressed by the Pauli matrices as 
(**,««) = \{cf",a z ). 

By the checkerboard decomposition |6| |j with the Suzuki- Trotter formula we express the 
partition function as 

Z = Tr[exp(-/3tf)] = lim Tr[cxp(-eiJ 1 )exp(-eiJ2)exp(-eiJ 3 )exp(-eiJ4)] m . (2) 

m — >oo 

Here (3 = 1/T is the inverse temperature, e = (3/m and the set of Hi arises from a decom- 
position of H as defined in ref. |l0[ . Since the coupling constant J can be absorbed into 
the temperature, we put J = 1 in our later discussions. Now we insert complete sets of 
eigenstates | + 1} and | — 1) of a z between each factor of exp(— eHi) to map the 2-dimensional 
quantum system to an induced (2+l)-dimensional classical system with Ising-like variables 
a(i,r) = ±1/2, 

Z= ]T exp[-S({ S (i,r)})]. (3) 
g(t,r)=±l/2 

Here r labels the slices in the artificial third dimension which has total 4m layers. S({s(i, r)}) 
consists of the four-spin interaction associated with the r-like plaquettes, which are called in- 
teracting plaquettes. Each interacting plaquette is bounded by four non-interacting plaquettes 
to form a checkerboard lattice |l2]| . The interacting plaquette weights are determined by the 
products of the transfer matrix 

exp[-S p {sx, s 2 ; s 3 , s 4 )] =< si, s 2 |exp[|(cr? ■ a * + a\ ■£7^)]|s 3 ,s 4 >, (4) 
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for each plaquette configuration C p (si, s^; S3, S4) with si and S2 locating on a rth-slice, 
and S3 and S4 on the (r + l)th-slicc [Q. Actually S({s(i, r)}) = X) P ^p( s ii s 2; S3, S4) and 
the sum is over all possible interacting plaquetts p. A calculation of the transfer matrix 
exp[—S p (si, S2; S3, S4)] leads to the form 



exp(|) 



/cosh(f) 

exp(— |) cosh(|) exp(— |) sinh(|) 

exp(— |)sinh(|) exp(— |) cosh(f) 

Vsinh(f) 




(5) 



The states of the spin-pair (si,S2) for the row index and (s3,S4) for the column index are 
arranged in a sequence of (+, +), (+, — ), (— , +), (— , — ). The elements of the matrix (up to 
an overall constant) can be interpreted as the Boltzmann weights of the configurations for 
the induced classical spin system. The zero elements indicate that only a part of the spin 
configurations is allowed. 

Updating schemes based on the transfer matrix must satisfy the ergodicity and detailed 
balance condition. Due to the facts that any allowed configurations can be achieved by flipping- 
pairs of spins on plaquette edges from an allowed configuration and each spin is shared by two 
interacting plaquettes it follows immediately that one should flip a closed loop of spins, i.e. 
change the signs of all spins on the loop simultaneously. Practically, we adopt the updating 
procedure introduced in refs. j^, [l3| in our simulations. A complete Monte Carlo sweep 
updates: 1) all spatial non-interacting plaquettes bounded by four interacting plaquettes; 
2) all r-like loops formed by three connected non-interacting plaquettes bounded by eight 
interacting plaquettes; 3) all Polyjakov lines in the r-dircction; 4) all global loops generated 
by connecting diagonal spin-pairs in the interacting plaquettes and vertical spin-pairs in the 
non-interacting plaquettes for all (i,r)-like planes. 

We first consider a dynamic relaxation process starting from an ordered state (mo = 1) 
[ p2[ |23| . In this dynamic process, it is believed that at the transition temperature Tkt or 
below, there exists a dynamic scaling form, e.g. for the kih moment of the magnetization, 

M< fc >(t) = b~ kri ' 2 M^\b-%b- l L), (6) 

which sets in right after a microscopic time scale t m i c . This time scale can depend on initial 
conditions, algorithms or other microscopic details. In simulations of classical spin systems 
tmic ~ 100 Monte Carlo time steps is observed. If a Monte Carlo time step is considered to be 
a typical microscopic time unit, this is reasonable. 

Starting from all spins up, we have performed the simulations with the Metropolis algorithm 
at the transition temperature Tkt — 0.3423 jy]]. Lattice sizes are taken to be L? x Am = 
64 2 x 120 and L 2 x 4m = 32 2 x 120. Samples for average are over 1 000 for L = 64 and 
2 000 for L = 32. Statistic errors are estimated by dividing the samples into three groups. We 
measure the magnetization defined as 

M(<) = Z^4^ < ^> (l ' r ' t)> (7) 

(i,r) 

and its second moment. From the scaling form in (^), it is easy to deduce that the magneti- 
zation decays by a power law 

M(t) - r" /2z . (8) 



Such a power law decay has actually been known for a long time in the long-time regime of the 
dynamic evolution, but now it is expected to hold also in the macroscopic short-time regime. To 
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determine the dynamic exponent z independently and further confirm the short-time dynamic 
scaling, we introduce a dynamic Binder cumulant U(t,L) = M^ 2 '/M 2 — 1. Simple finite size 
scaling analysis shows 

U(t,L)~t d/z . (9) 

This behaviour of the Binder cumulant is a typical behaviour in the short-time regime, where 
the non-equilibrium spatial correlation length is very small. 
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Fig. 1. - The time evolution of the magnetization for L — 64 and L = 
ordered state is plotted in double-log scale. 



32 by starting from 



Fig. 2. - The time evolution of the cumulant for L — 64 and L = 32 by starting from ordered 
state is plotted in double-log scale. 



In fig. 1 the time evolution of M(t) is plotted in double-log scale. The solid line is that for 
L = 64 and the circles for L = 32. One can see that until 1000 Monte Carlo steps, there is no 
visible finite size effect for these lattice sizes. Here the microscopic time scale t m i C ~ 100. In 
the figure we have skipped data for t smaller than 200. From the slope of M(t) in the time 
interval [200, 1000] we measure the critical exponent rj/2z = 0.0522(6) for L = 64 and for 
L = 32 we get rj/2z = 0.0518(7). Within the statistic errors, they coincide. 

From the measured M(t) and M^ 2 \t) we calculate the Binder cumulant U(t) which is 
plotted in fig. 2. It is interesting that the curves of U(t) exhibit short-wave fluctuations, which 
have not been observed in classical spin systems [^2| p3| . However, the short-wave fluctuations 
do not affect so much the long-wave behaviour and a clear power law behaviour is seen in 
the figure. From the slopes of the curves, we measure the critical exponent d/z = 0.85(3) 
for L = 64 and d/z — 0.86(3) for L = 32. With the results of L = 64, we get the critical 
exponent z — 2.35(8). Taking z as input, from the measured value of r]/2z — 0.0522(6) we 
obtain r\ = 0.245(8). Compared with the results from simulations in equilibrium, our value of 
the exponent t] supports r\ — 0.25 ± 0.01 in ref. Q but is smaller than r\ — 0.29 ± 0.09 or 
r\ = 0.276 ± 0.014 given in ref. ||. This suggests that the exponent rj in this quantum XY 
model takes a classical value r\ = 1/4. 
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In the short-time regime of the dynamic evolution, there are plenty of new phenomena. 
To describe all these phenomena, the static exponents together with the dynamic exponent 
z arc in general not sufficient. Another interesting and important dynamic process is the 
relaxation starting from a high temperature initial state with a small initial magnetization 
mo. It has been shown analytically and observed numerically in classical spin systems, that 
after a microscopic time scale t m i c , the magnetization undergoes surprisingly a power law 
initial increase JIJ, [lj| [l7| 

M(t)~m Q t e . (10) 

The exponent 9 is a new independent critical exponent, which can not be expressed in terms 
of other known exponents. 

In order to measure the critical exponent 9, we must prepare an initial state of random 
configurations but with a small initial magnetization tuq. In the high temperature limit of 
the transfer matrix in eq. (^), the four elements proportional to sinh(|) vanish and only the 
four equal-weighted non-zero diagonal elements remain, which correspond to the interacting 
plaquettes of C p (si, s 2 ; s 3 , s 4 ) = (+,+;+,+), (-,-;-,-), (+,-;+,-) and (-,+;-,+). 
Therefore, we first put randomly s* = h with a probability of p = (m + l)/2 and s\ = — h 
with a probability of 1 — p respectively on the first layer (r = 1) of the L 2 x 4m lattice. Then 
we copy this configuration to all other r — 1 layers. This procedure gives the required initial 
random configuration. 

The simulations have been performed with a lattice size L — 64 and total samples over 
4 000. It is well known that the microscopic time scale t m i C in this case is small, typically less 
than 20 Monte Carlo time steps |l6|, [l7| [2l], ^3). Therefore we stop updating at t = 150. In fig. 
3, the time evolutions of M(t) are plotted in log-log scale. After about 30 Monte Carlo steps, 
the curves show a nice power law behaviour. From the slope of the curves, we measure the 
critical exponent 9 = 0.18(1) for mo = 0.01 and 0.18(2) for mo = 0.005. Rigorously speaking, 
the exponent 9 is defined in the limit mo — > 0. However, within statistical errors our results 
for mo = 0.01 and 0.005 show already no difference. Therefore an extrapolation of 9 to the 
limit mo = is not necessary here. In conclusions, for the first time we have investigated the 
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Fig. 3. - The time evolution of the magnetization for L = 64 with mo = 0.01 and mo = 0.005 
is plotted in double-log scale. 



short-time critical dynamics of the quantum spin system taking the two-dimensional spin 1/2 
quantun XY model as an example. Universal power law behaviour of the magnetization and 
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the Binder cumulant is found and the new critical exponent 9, the dynamic exponent z and the 
static exponent r\ are determined at the transition temperature Tkt — 0.3423. The measured 
exponent r\ = 0.245(8) agrees very well with 77 = 0.244(5) obtained for the classical XY model 
[ p2| and also with the theoretical prediction 77 = 0.25. However, the exponent 6 = 0.18(2) and 
the dynamic exponent z = 2.35(8) are apparently different from 9 = 0.250(1) and z = 1.96(4) 
for the classical XY model jnl ^2| . These results indicate that the quantum XY model and the 
classical XY model are in a same static universality class but in a different dynamic universality 
class. As an alternative approach, our simulations pave a way for a systematic application of 
the short-time dynamic scaling to quantum spin systems. 

For quantum spin systems, besides the stochastic relaxational dynamics there exists a true 
dynamics described by the Heisenberg equation of motion. It is challenging whether the 
short-time dynamic scaling can be found also there. 
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